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ABSTRACT 


We investigated an inverse imaging problem by applying the Rytov 
approximation for Through-the-Wall Imaging using THz waves. In the beginning, we 
studied some properties of THz waves and the physical conditions for THz imaging in 
matter. Then we started with Maxwell’s equations to derive a model for the transmission 
of Green’s functions and used a Lippman-Schwinger integral equation and Rytov 
approximation to predict the scattered field. We applied the L-curve method for the 
selection of regularization parameters, and then presented the reconstruction algorithm 
and illustrated the result with numerical simulations. We also compared this result to the 
one obtained by the Bom approximation. 
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I. INTRODUCTION 


The amazing ability of terahertz radiation to penetrate through materials allows 
for the nondestructive, noninvasive inspection of mail and packages in post offices, and 
luggage at airport security check points [1]. 

The term terahertz (THz) refers to radiation whose spectrum lies between the 
infrared and microwave bands with frequencies ranging from 0.1 to 10 THz, where 
lTHz= lO'^ Hz or, in units of wavelength, X =3 mm to 30 //m. Until recently, this “THz 
gap” has been almost inaccessible due to the lack of efficient sources and detectors in this 
region [2]. But recently, sources and detectors for THz frequencies have become 
available and easier to use. Therefore, many applications of terahertz radiation, such as 
imaging and communications, are now appearing. Among these applications, a very 
important application (and the focus of much research) is imaging using THz radiation 
[3]. The study of this imaging is an interesting and attractive area for researchers. 

A. BACKGROUND 

Various THz imaging techniques have been developed since Hu and Nuss first 
demonstrated THz imaging in 1995 [4]. THz imaging applications are the focus of 
constantly growing interest. Because THz waves can provide great spatial resolution and 
can penetrate most dielectric materials, such as plastics, paper, and wood, there are many 
applications of THz imaging. A classical application of THz wave based remote sensing 
is content inspection of packages and envelopes. Recently, THz imaging has been used 
for medical and biological applications because it provides high quality images. Another 
application is real-time nondestructive testing, which includes testing for defects in 
plastic food and medicine packages. Also, there are a number of security applications for 
THz imaging: for example, Through-the-Wall Imaging (TWI), luggage inspection, and 
the detection, from a distance, of weapons or explosives hidden under clothing or 
briefcases. 
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B, MOTIVATION 

In the decade between 1996 and 2005, more than a half million servicemen were 
injured or killed in situations involving barricaded offenders, hostage-taking, and high- 
risk entry [5]. Through-the-Wall Imaging (TWI) has been studied for saving soldier’s and 
officer’s lives. TWI technology can help soldiers to detect enemies hiding in a building, 
and, for this reason, TWI technology needs to evolve and improve. 

The armed forces personnel or police officers who experience operations on 
urbanized terrain may want equipment that can detect, identify, and localize enemies 
through a wall, and operate day and night in all weather conditions. 

Most imaging systems use bi-static or multi-static transmission, where the 
transmitter is on one end, and the receiver is on its opposite end. Mono-static systems, on 
the other hand, use transmitters and receivers that are in the same location. Mono-static 
imaging systems rely on the back-scattered field. Through-the-wall imaging is more 
difficult and challenging than any other back scattering method because of the general 
variability of wall construction. In this thesis, we will concentrate on mono-static imaging 
techniques. 

There are two basic problems in ray-scattering. The first one is the direct 
scattering problem, which predicts how the rays scatter from a known object. The second 
is the inverse scattering problem, which attempts to recover the original object from the 
blurred and corrupted measurements of the scattered field. The inverse scattering problem 
is one of determining the characteristics of an object from measurement data of radiation 
or particles scattered from the object. Also, most inverse problems are ill-posed. X-ray, 
Tomography, and ultra sound imaging are good examples of ill-posed inverse problems. 

Let’s look at the difference between the direct problem and the inverse problem. 
The forward problem can be expressed as 

D = A{f) (1.1) 


in which, we denote the data by D e □ “ , the object by / e □ ^ and the operator acting 
between two spaces containing D and f hy A. {N and M denote the dimensions of, 
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respectively, the object and data spaces and, in general, N ^ M .) So the direct problem 
predicts the data (D) from the object (/). On the other hand, the inverse problem is 
expressed as 

f=A-\D) (1.2) 

and solves for the object (/) from the data {D). This equation looks easy to solve, but 
real-world inverse problems are typically much more complicated. 

Now let’s look at the inverse problem a little bit deeper. If we want to accurately 
map objects to data, we must also consider the measurement errors. The difference 
between the actual measurement data (D) and the real data (D) is denoted by noise 
n = D-D, and the scattering problem can be more accurately expressed by 

D = A(f) + n (1.3) 

The object and data functions are conveniently considered as belonging to Hilbert 
spaces. In equation (1.3), A is an operator from one Hilbert space to the other. Hilbert 
space allows simple geometric concepts like projection onto fewer dimensional spaces 
and loss of information. 

Inverse problems are almost always ill-posed problems. J. Hadamard formalized 
the concept of well-posedness for equations of the form (1.3) [6]. Any equation that does 
not satisfy all three conditions is called ill-posed, and a method for solving this problem 
approximately is called a regularization method. 


1 

Existence 

There exists a solution to the problem 

2 

Uniqueness 

There is at most one solution to the problem 

3 

Stability 

The solution depends continuously on the data 


Table 1. Conditions for Well-posed Problems. 
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c. 


PROPERTIES OF TERAHERTZ FREQUENCIES 


The term “Terahertz” is applied to broadband pulsed electromagnetic radiation 
with a spectrum that falls between the infrared and microwave bands ranging from 0.1 to 
10 THz [7]. The wavelength in this region ranges from 3 mm to 30 jum (as shown in 
figure 1). These are short wavelengths, and they have high resolution compared to other 
radio frequencies [2]. There are several distinct advantages of THz waves. The first is 
that most chemical substances possess characteristic absorption features in the THz range. 
This property enables us to detect illegal chemicals, even when they are sealed inside a 
packet or clothing. The second advantage is material transmission properties. Most non¬ 
conducting materials such as plastics, paper, wood, and ceramics are totally or partially 
transparent in the THz range. The third advantage is that the radiation is non-ionizing and 
has relatively low energies compared to X-rays; therefore, we don’t have to be concerned 
about the safety of THz imaging. 


In this thesis, we will use THz pulses for imaging. Figure 2(a) shows examples of 
THz pulses transmitted through nylon. When a THz pulse interacts with an object, it 
experiences a delay, an attenuation, and a broadening, as the different component 
frequencies are phase shifted, absorbed, reflected, and scattered. Figure 2(b) shows the 
power spectrum that the pulse frequency content is centered around ITHz (note that the 
higher frequency content is attenuated more than the lower frequencies by nylon). 



Figure 1. Electromagnetic spectrum[From 7] 
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Figure 2. Examples of terahertz pulses in (a) Time Domain and (b) Magnitude of 

Fourier Coefficients [From 4] 

For the convenience of modeling, we assume that the wall that we want to see 
through is a linear, homogeneous, and isotropic dielectric. We know THz frequencies 
have good penetration for dry materials, but a spectral cut-off above 3 THz will be caused 
by absorption from water vapor in the THz beam path [8]. So we are not generally 
interested in frequencies over 3 THz. We also need to consider the object that we want to 
image. The reflection from the object depends on the object properties. For example, 
metallic materials reflect perfectly, but skin or water absorbs THz rays, so we can’t get 
significant reflected pulses. Therefore, we should be careful to choose frequencies that 
balance between the penetration of medium and the reflection of the object. 

D. THESIS ORGANIZATION 

This thesis uses Born and Rytov approximations to model measured scattered 
fields from objects located behind walls. 
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Chapter II introduces the concept of inverse problems, including ill-posed 
problems, singular value decomposition, and regularization methods. The Born and 
Rytov approximations are discussed as well. 

Chapter III provides some properties of electromagnetic waves in media and 
restricts our model system to a real situation. A scalar theory is introduced for use in the 
succeeding chapters. 

Chapter IV develops the transfer function of the media. This chapter describes the 
mathematical approach and derives Green’s functions in different regions. 

Chapter V examines reflectives of media and evanescent waves. It also examines 
materials with reflectivity. 

Chapter VI derives a wave equation in the specified medium and the scattered 
fields by applying the Born and Rytov approximations. 

Chapter VII provides algorithms and simulations (with noise) by applying 
Tikhonov regularization and Truncated Singular Value Decomposition regularization 
schemes. 

Finally, Chapter VIII gives a final summary conclusion, comparing the simulation 
results and providing suggestions for future work. 
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II. THE INVERSE PROBLEM 


This chapter will introduce ill-posed problems and regularization solution 
methods more deeply. Regularization methods are used to mitigate noise effects in the 
data or the operator. Three widely used regularization techniques are Tikhonov 
regularization, Truncated Singular Value Decomposition, and regularizing iterative 
methods. In this thesis, we will use the Tikhonov regularization method and the 
Truncated SVD method. We will confine ourselves to linear equations with compact 
operators in Hilbert space and use the Singular Value Decomposition for our 
reconstruction algorithm design. 

A, REVISIT ILL-POSED PROBLEMS 

We want to get images using electromagnetic frequencies. So there is no way to 
avoid the inverse problem. The very first equation we need is a Lippman-Schwinger 
equation [9]. 

We denote: 

• G^(x,x') is Green’s function which satisfies the wave equation with an 
impulsive source term 

• p(x') is a source factor determined by the index of refraction of the 
scatterer 

• i//(x) is an electric field, tPi„^(x) is the incident field, and is 

the scattered field 

Equation (2.1) is derived from the reduced scalar wave equation 

V^y/{x) + en\x)y/{x) = Q (2.2) 

where n(x) is the index of refraction and k is the wave number. 
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In this thesis, we ignore multiple scattering. Because we assume multiple 
scattering terms are very small in comparison to primary scattering, terms we can rewrite 
equation (2.1) as 

¥sca,A^) = (2.3) 

To get images of objects, we have to estimate yC>(v) (which defines the “image”), 
but there are some problems. The main problem is that we can never obtain data for all 
points (v) and at all frequencies {co). The other problem is that the measurements are 
always accompanied by noise. If we consider these problems, we can make a more 
complete measurement model and express equation (2.3) in matrix form. To perform this 
discritization, we break the integral in equation (2.3) interval x e [a b\ into M parts and 
X e [c d] into N parts. Then the integral equation becomes 

fb 

¥s.{X:S)f{s)ds = d{x.) a<s<b, c < x. < X 2 < ■■■ < Xj^ < d (2.4) 

J a ’ 

where operator K is defined as the kernel which operates on the object f e X to 
produce the measurement d . We have simplified equation (2.3) into a one 
dimensional case and changed notations to better understand this measurement system. 

Where, 

• The measurement at position X; : ^i^^_.^„(X;) ^ d(x.) 

• The “kernel” : (x,x') ^ K(x,.,5) 

• The “object” function : p(x') ^ f(s) 

Then we can rewrite equation (2.4) in matrix form (with noise included)as 

d = Kf + n (2.5) 

where 
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J(Vj) 


K(xi a -1- As) 

K(xj a -1- 2As) 

■■ K(xia + MAs) 

d = 

dix^) 

K = 

Y^ix^a^- As) 

K(v 2 a + 2As) ■ 

■■ Y^(x2a+MAs) 


_d(Xj^)_ 


K(v^ a -1- As) 

K(v^ a -1- 2As) ■ 

■■ Y^(Xj^a + MAs) 



f(a + As) 


n(Xi) 

f = 

f(a + 2As) 

, and n = 

^(Vj) 


f(a + MAs) 


_n(v^)_ 


and d e and / e are vectors. 

Let’s take a close look at the dimension of the measurement space and the object 
space. The measurement space has dimension N , and the object space has dimension M . 
This reflects the number N of the discrete measurements and the number N defining the 
resolution with which / can be estimated. Normally the object dimension M is greater 
than the measurement of dimension N. Now, equation (2.5) can be used to get / by 
defining an inverse K ^. However, the noise can alter the measurement space and lead to 
an unstable estimate of / = K ' J . Equation (2.5) is also often an over-determined system, 
which means there are more equations than unknowns, and we can estimate / by least- 
squares. Therefore, we have to figure out an over-determined system to get solutions 
close to the real (discrete) values. We will deal with some methods to get stable solutions 
in section D. 


B, LEAST SQUARES SOLUTION BY SINGULAR VALUE 
DECOMPOSITION 


We consider d to be a vector with components d^ =d(x-). Normally, the vector 

of measurements d in equation (2.5) will not lie exactly in the measurement space 
Y (because of noise). Because J e T only if the noise vectors e Y and Kf eY. But n is 
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Figure 3. Hilbert space and measurement space 

independent of the operator K , and the object vector / may have components in 
measurement space or in null space. So a possible “best object” / will minimize the 
distance \d -I^|| • We denote the least squares solution to equation as 

/ =arg min||j-K/|| (2.7) 

From Figure 3 we can get the solution to equation (2.7). The vector K/ lies in 
measurement space Y . When ¥^f = P , Kf gets the closest distance to d . The vector 
d - P = d - Kf is orthogonal to the measurement space Y and all vectors of the form Ky 
in measurement space Y. Then we can get the normal equations by [10] 

K^Kf = K'^d (2.8) 

/ is the least square solution to equation (2.8). Since the matrix K^K e □ is square 
we can write 

7 = (K^K)-'K^d (2.9) 

/ is exactly what we want. In order to get / , we have to compute (K^K) . But this 

can be painful to invert with large N and M . So we want to use another method, 
singular value decomposition, to get / . 
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We would like to employ simple linear algebra to build matrices of the form (K^K) . 

Consider a matrix A that is square (MxM) with an eigenvector, Xj an eigenvalue. 

We can transform the matrix A in terms of the matrices e-e[, formed from the set of 
vectors {ej, using the Spectral Theorem which states that A can be written with as 

M 

A = UAU-^=UAU^=Y,\e,e] ( 2 . 10 ) 

1=1 

where U is the matrix with the eigenvectors of A on its columns and A is the matrix 
with the eigenvalues of A along its diagonal. 

We assume is an orthonomal basis of □ “ , and so we can represent an 
M dimensional vector y in terms of e. 

M M 

y = ^( 4 yX=^yA ( 2 . 11 ) 

/=i (=1 

where y. = ef y is the component of the vector in the basis. Consequently, 

M MM 

X = Ay = ^ y = y,e, = £ x,e, ( 2 . 12 ) 

/=1 (=1 /=1 


This is representation of x = Ay in terms of the basis {ej. 


Let’s adjust this result for the matrix(K^K) . Suppose are the eigenvectors of 

K^K. Then 


(K^K)^‘ =([/At/^')^' 


( 2 . 13 ) 


Multiply by on the right of both sides 
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(K^Kr'K^ =t/A^'(K[/f 



r 




UA^‘ 

1 

Kgj Kcj • 

1 

... 

y 


1 

> 

0 

1- 

0 

-1 


o • • 

(J 2 ! \ 

0 • • 

... 


1 

o • 

■ 0 

_ 

... 



D yT 


= UDV^ 


We define cr,. = Ke^ and cr,. = 


jo-r^Ke,. i/cr,. ^0 

[O ifcT,=0 


Now multiply both sides of equation (2.9) by K so that 

k 7 = K(K^K) ‘K^d = P 


(2.14) 


(2.15) 


P is the projection of the measurement vector onto measurement space. Therefore, 
K(K^K) is a projection matrix mapping d to its components in measurement space. 
We will use this projection matrix to get more accurate data using regularization 
techniques in the following section. 


C. TRUNCATED SINGULAR VALUE DECOMPOSITION 


To obtain a better estimate of the least square solution, the truncated singular 
value decomposition solution is often used. This method truncates the SVD 
representation to reduce the effect of noise contamination. Just select one small threshold 
value fC and compare this with singular values /I, .If /C > /I,., then we replace the l/^^ 

in D by 0. In this thesis, using the plot of singular values, we will find the critical point 
where the singular values change steeply and will truncate at that point [11] (see Figure 
21 ). 
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D, REGULARIZATION METHODS 


Since our model is linear and ill-posed, the inverse image is under-determined, 
and we will get small singular values of K. The problem is that K depends on our model 
of the measurement process and is not completely known. So we get singular values with 
a slight imprecision. Now we introduce regularization methods to overcome the problems 
associated with the small singular values that lead to reconstruction instability. 
Regularization is important for solving inverse problems, because the output of least 
squares is affected by data and rounding errors. Therefore we would like to introduce 
regularization methods to reduce these errors. But there is a trade-off between the 
regularized output and the original sets of data. If the regularization is too crude, the 
regularized solution does not fit the given data, and the residual error ||r|| = ||d-I^|| is 

unnecessarily large. If the regularization is too small, the fit will be good, but data noise 
effects will be larger. 

There are also some regularization methods that employ operator or data 
correction [12]. 

1, Regularization by Operator Correction 

Let’s model the problem by mapping A as 

Ax = d, xeX, deY (2.16) 

The equation (2.16) is said to be an ill-posed problem, if it does not meet one of the 
conditions in table I. If we want to transform an ill-posed problem into a well-posed 
problem by approximation of the equation, then we have to choose a mapping method. 

The regularization of equation (2.16) consists of a substitution of the operator A 
by A . Let A • A ^ T be an operator such that 

Ax = d (2.17) 

This is the idea behind the well known Tikhonov regularization. 
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2, Regularization by Data Correction 


Let K Y ^ A be a continuous operator such that 

then 

Ax = Kd, x^X, d^Y (2.18) 

is well-posed and called as regularization by data correction. 

3, Regularization by Data and Operator Correction 

The last method of regularizations is a combination of both previous methods. Let 
A • Z ^ T be an operator and K Y -^Y a continuous mapping such that 

Ax = Kd, x^X, d^Y (2.19) 

Tikhonov regularization is also a well known method of this kind of regularization. 

E, L-CURVE METHOD 

Since the quality of the regularized solution is controlled by the regularization 
parameter for example, the threshold parameter /C used in truncated SVD regularization 
[13], we have to choose the optimal parameter to get a “best” image. There are various 
methods for the selection of regularization parameters, including the Discrepancy 
Principle, Generalized Cross Validation, and L-curve. The L-curve method applies a log- 
log plot of the regularized solution against the squared norm of the regularized residual 
for a range of values of regularization parameters [13] and selects the corner as the point 
of maximum curvature in the L-curve plot (see Figure 4). 

There are two main approaches used in the L-curve method. The first approach 
considers both the residual norm and the solution norm. The second approach considers 
the value of the maximum curvature. Like other methods, the L-curve method has its 
merits and limitations. The merits of L-curve are that it’s robust and can treat 
perturbation caused by correlation noise. On the other hand, the L-curve method is 
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limited, because it is not convergent [14], i.e., since singular values are getting close to 0 
but not exactly 0 as increasing index of singular values, we can’t get the value which 
satisfies the norm is 0. 



Figure 4. L-curve for Tikhonov Regularization [From 14] 


F. APPROXIMATIONS TO THE WAVE EQUATION 

In the case of electromagnetic scattering problems, we need to approximate the 
unknown internal field by a known distribution in order to linearize the problem. The 
Born and Rytov approximations are widely used to simplify electromagnetic scattering 
problems. In this thesis we will derive reconstruction formulas under the Rytov 
approximation and compare with the Bom approximation. 


I, Born Approximation 


The first Born approximation was introduced by Max Born in 1925, in order to 
solve problems concerning the scattering of atomic particles [7]. The Born approximation 
is obtained by expressing the total wave field as the sum of incident field plus scattered 
field: 


y/{r)=y/^{r) + y/Xr) 


( 2 . 20 ) 
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This approximation assumes that scattered fields inside scatterers are negligible 
compared to the background and the normal fields. The first-order Bom approximation is 
good only if the scattered field is much smaller than the incident field. 

2, Rytov Approximation 

The Rytov approximation was discovered in 1937 by Rytov, who was analyzing 
the diffraction of light by sound waves. This approximation was later applied by 
Obukhov to describe the propagation of electromagnetic waves in random media [15]. 
The Rytov approximation is based on expressing the total field as a complex phase, 
which is the sum of incident phases plus scattered phases: 

= ( 2 . 21 ) 

The condition for the applicability of this approximation is that the phase of the scattered 
field varies slowly, relative to one wavelength. So the size of target is less critical in the 
Rytov approximation. 
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III. ELECTROMAGNETIC WAVES IN MEDIA 


It is difficult to explain electric fields and magnetic fields, which are cross- 
coupled, without the Maxwell’s equations. Maxwell’s equations were established by 
James Clerk Maxwell by 1864 and experimentally verified by Heinrich Hertz in 1888 
[16]. Since then these equations have been applied in a variety of areas such as optics, 
microwaves, antennas, radar, and communications. In this chapter we will investigate 
Maxwell’s equations in media and derive the associated wave equations. 


A. MAXWELL’S EQUATIONS 


Maxwell’s equations are comprised of Gauss’ law for the electric field, the 
observation that there are no magnetic monopoles. Ampere’s law, including the 
dD 

displacement current —, and Faraday’s law of induction in differential form. They can 

dt 

be written as 


VD = p 
V-5 = 0 


VxE 


dB 

dt 


VxH=J+ 


dP 

dt 


(3.1) 


where D and B are the electric and magnetic flux intensities, and E and H are the 
electric and magnetic field intensities, J is electric current density, and p is volume 
charge density. In reality, there are two more things we have to consider for Maxwell’s 
equations. These are P , and M , called induced polarization and magnetization 
respectively. The associated constitutive relations are expressed as 


T) — ^()E “1“ P 

(3.2) 

5 = //„(//+M) 

(3.3) 
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Using equations (3.2) and (3.3), we can rewrite Maxwell’s equations in terms of E and5 . 


V ■E = — {p-V-P) 

% 

VB = 0 


VxE = 


dB 

dt 


BE 

V X 5 = //()Gq + //,) 

dt 


dP 

J + — + VxM 
dt 


(3.4) 


where and e,, are the permeability and permittivity in free space, respectively. The 
speed of light (c) and characteristic impedance (; 7 „) of a vacuum are expressed using ju^ 
and e„ as 


c = 




(3.5) 


1, Simplified Model 


Electromagnetic fields behave differently in the presence of media. For Through 
the Wall Imaging using electromagnetic waves, we will confine our attention to materials 
obeying a linear model. Materials used for constructing walls include wood, concrete, 
and iron-beams. These materials typically have anisotropic permittivity. Anisotropic is an 
inherent property of the atomic and molecular structure of the dielectric [17]. So in 
anisotropic materials, permittivity depends on its direction and the linear constitutive 
relation can be written in matrix form as 



(3.6) 


In general, the field vector E is no longer parallel to D . But we can simplify our model 
by using the principal axis system: We consider homogeneous and linear materials so that 
the permittivity is independent of position. Then the permittivity tensor can be expressed 
with six elements as 
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^xy ^xz 

e e e (3.7) 

e e e 

xz yz zz _ 

Since this matrix is symmetric, its eigenvectors are orthogonal. If we choose a 
coordinate system which is aligned to these eigenvectors, we can rewrite the permittivity 
tensor as a diagonal tensor determined by its eigenvalues. Then equation (3.7) becomes 

■e„ 0 O' 

tensor = 0 0 (3.8) 

_0 0 

Let’s investigate solutions to the Maxwell’s equations for source free 
homogeneous isotropic media. Here, J = p = 0 and Maxwell’s equations (3.4) become 


tensor(e) = 



(3.13) 



(3.14) 


v'£-//e^ = 0 

dt 

The plane wave solutions to equation (3.14) can be written as 

= (3.15) 


Substituting equation (3.15) into equation (3.14) we find the following equation 

(3.16) 


This is called as the dispersion relation. We can easily find the phase velocity as 

colk = ll4Jie (3.17) 

For conductive materials with permittivity e+ icTjco, the wave number is complex and 
written as [16] 


k — kjf + ikj 


~ CO 




cr 


coe 


1/2 


= ( 1 + 0 , 


cojua 


(3.18) 


Then the penetration depth is 

2 

cojua 

For a wave propagating in the +z direction, we have z dependence given by 

^-ikz _ ^-ikgZ-ik,z 



(3.19) 


(3.20) 


This wave attenuates exponentially in the same direction that it propagates. Figure 
5 shows a wave propagating in a lossy medium. 

1 cr 

l< — 

For a highly conducting medium with the penetration depth is a very 

small. So the waves will not be able to penetrate conductive materials. Therefore we 
ignore highly conductive materials in this thesis. 
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Figure 5. Wave propagating in a lossy medium[From 18] 

2, Properties of Frequency and Dispersion 

When the electric field is applied to a dielectric, the electric field tends to separate 
the electron from the positively charged nucleus, and this creates an electric dipole 
moment. Moreover, the dipole moment contributes to the electric flux density so that 

D = e„F + P = e„(l + = e{a))E (3.21) 

where j(o) is susceptibility. From equation (3.21) we can see the permittivity is 
changed because of the dipole moment. The polarization of a medium with refractive 
index n is 

P = X%E 

and the susceptibility is expressed as [16] 

2 

nq 

X ~ - - - 

me^icol -o/ -iy) 


(3.22) 


(3.23) 
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where / = j/ , which is a measure of the rate of collisions per unit time. Rationalizing 
Equation (3.23) as 


;r = 



(3.24) 


and substituting equation (3.24) into equation (3.21), we get the permittivity as 


^ inq^co/ 

o o\2 ooH F/ o o\2 


( 1 2\ 

^ , 2 2 

- + — 

(2 2\ 

^ , 2 2 

(®o ) 

+ CO / 

m 


+ CO / 


(3.25) 


When cl> is not close to , the imaginary term will be much smaller than the real term. 

This means e is approximately real valued when co is not very near resonance. When not 
close to a resonance, the deviation of the refractive index from unity is 


n-lU 


nq 


2me|) - 0 )^^ 


1 + 


ICO/ 


2 2 
COq—O) j 


(3.26) 


From equation (3.26) the real part of n represents dispersion and the imaginary part 
represents absorption. Figure 6 shows us that around the resonant frequency co ^^, the real 
part of n behaves in a strange manner and drops rapidly with frequency, and the material 
absorption occurs quite high. So we have to avoid choosing a radiation frequency that is 
in resonance with the dielectric. 


To make our model system more numerically tractable, we need more 
assumptions. Dispersive media have values of and ^ that depend on frequency. As 
a result, the wave velocity is generally frequency-dependent [18]. So, dispersive media 
have frequency-dependent constitutive laws; this means waves of different frequencies 
propagate with different velocities. For dispersive media, we need lots of information 
about different materials. The analysis is complex and unnecessary. So we will not 
consider highly dispersive media in this thesis. 
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Figure 6. 


Real and imaginary parts of the refractive index [From 19] 


In general, the velocity of the envelope of any modulated sinusoid is the group 


velocity ^ which is the speed at which energy and information travel [18]. 


v„ = 


\da) j 


(3.27) 


For these reasons, our model assumes that the media have low conductivity and are non- 
dispersive. 


B. THE SCALAR THEORY 

In the previous section, we derived the vector wave equation (3.14). Also, we 
restricted our model to a dielectric medium that is linear, homogeneous, isotropic, and 
non-dispersive. If the medium is just like that, we can rewrite the wave equation for the 
electric field as 

V"E-^^ = 0 (3.28) 

c dt 

Identically, we get the same equation for the magnetic field as 

V'H-^^ = 0 (3.29) 

c dt 

Since these are vector equations, we will get 6 equations with 6 components: E^,E^,E^ 
and H ^,H y,H But fortunately, we can make these equations into a single scalar wave 
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equation because under our model restriction, all vector components of E and H are 
identical. A single wave equation becomes 


rv 2 / ^ d^y/{x,t) _ 

vVU0+^ A -0 


(3.30) 


where represents any of the vector field components and y/ is dependent on 

position X and time t. However, when the medium is inhomogeneous with a permittivity 
e(v) that depends on position v, then the wave equation is not the same as before. The 
wave equation is derived by the following procedures. First we consider the first line of 
equation (3.1) and equation (3.21). Combining these equations yields 


eV-E + (Ve)-E = p 


(3.31) 


Divide by e on both sides 


V E^P_(Zl:^^£_(Vine)-E 

e e e 


(3.32) 


Now take curl of the third of equation (3.1) 


Vx(VxE) = -Vx — 
dt 


(3.33) 


V(V-E)-V"E = - —(Vx(//H)) 


(3.34) 


And substitute equation (3.32) into equation (3.34) 


V ^-(Vlne)-E -V"E = -—(Vx(//H)) 


(3.35) 


One of the laws for operations with the div-and-curl operator is 


V x{(pF) = (f>V X F + (V (f>)x F 


(3.36) 


and so equation (3.35) becomes 


V"E-V^ + V[(Vlne)-E]= // —(VxH) + —[(V//)xH] 
© dt dt 


(3.37) 
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Using the fourth line in equation (3.1), equation (3.37) becomes 

V'E-Qu^ + [V(Vlne)-El = V^ + —[(V//)xHl (3.38) 

dr ^ ^ e 5r ■' 

For nonmagnetic material, ju = ju^. When p = 0 , equation (3.38) can be simplified to 

02-17 

V"E + V[(Vlne)-E]-e//^ = 0 (3.39) 

Finally, with the help of the refractive index, = ^ , and the speed of light in vacuum, 

/ ®() 

^ , equation (3.39) becomes 

V"E + 2V[V(lnn)-E]--^-^ = 0 (3.40) 

The extra term in equation (3.40), 2V [V(ln n) • E] , will not be zero when n varies with 

position. So electric fields and E^ may not satisfy the same wave equation. 

Consequently, we can’t make these vector waves into a single scalar wave. Also, we will 
get some error by using scalar theory, even in the homogeneous medium when the 
boundary conditions are imposed on a wave [20]. 

We will consider scalar components and the corresponding scalar wave equation 
instead of the full vector theory in the thesis. 
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IV. THE TRANSFER FUNCTION 


For the purpose of our model, we have to develop the transfer function of a wall¬ 
like medium. We assume that the wave propagates into the medium with an index of 
refraction {n^). The source is located at = -ZqZ, and coordinates are as shown in 
Figure 7. 

The scalar wave equation with source s(r,t) in the time domain becomes 

V^E-4^ = s(r,t) (4.1) 

c dr 

Using the Inverse Fourier Transform, we can express the source term in terms of angular 
frequency ( co) and oscillation frequency (v) as 

1 /•OO 

s(r,t) = — [ Sir,(D)e'‘^‘dco (4.2) 

2;r 

where S(r, co) is the source term in the angular frequency domain. 

Using the identity co = Inv and da) = Indv , equation (4.2) becomes 

s(r,t)=r S(r,v)e‘^^^‘dv (4.3) 

J -00 

where S(r,v) is the source term in the oscillation frequency domain. 

In the previous chapter, we showed that a plane wave solution is in the form of equation 
(3.15). This will satisfy equation (4.1), and we can obtain the following relationship. 



Then we can make equation (4.1) into a time independent wave equation as 

V^E + n^eE = S(r,v) (4.5) 

To get the solution of equation (4.5), we put the source at the point = -ZqZ and use the 
Green’s function which satisfies the following equation 
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Figure 7. Point source emanating waves through a medium 


(V^ + k^n{z))g{r,r^) = 5{r + r^) = 5{x)5{y)5{z + Zq) 


(4.6) 


where k = and /I is the wavelength in free-space. 


The refractive index of our system is 


n{z) = 


1 for z < 0 
n for 0< z< L 
1 for L< z 


Then we can get the solution for E in the frequency domain as 

£(r,v) = j g{r,rf)S{r,v)d^r^ 


(4.7) 


(4.8) 


Where g{r,rf) is the Green’s function in the spatial domain, and the solution for E in 
the time domain (using the Inverse Fourier Transform) as 


E{r,t) = r E{r,v)E^^''‘dv 

J -00 


(4.9) 
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In our model, we assume that the medium is essentially non-conductive, so the 
electromagnetic wave can propagate the medium with little damping. Since we put the 

source at “ ° along the z direction we can convert the equation (4.6) into the 

Fourier domain. 


/•OO /“OO 

G(v^,v^,z,r^,)=\ gix,y,z,r^,)e^ 

J -00 J -00 


(Ittv^x+Ittv y) 


dxdy 


(4.10) 


gix, y, z,r^,) = ^ ^G(v^,v^,z, dv^dv^ 


(4.11) 


These expressions are correct in our model, because the medium remains constant in the 
X and y direction. Therefore, there is no changing of the medium in these directions. 


Now substituting equation (4.11) into equation (4.6) and using the identity in 
equation (4.10), we obtain 








G(y,,v, 


z,r^^) + k^rdiz)Giv^,v^„z,r^^) = ^iz + ZQ) 


(4.12) 


If we define 

F = 2nv^x + 2nVyy 

we get 

\Ff ={27cvf + {2nv^f 

Using the wave number definition, we can rewrite 


e=ki+ki+e = 

X y z 


^2;r^ 

2 


2 

'^2;r'' 

J 

+ 


+ 

iXj 


= ( 27 cvf +( 2 nvS +( 2 nv? 




(4.13) 


(4.14) 


(4.15) 


If z^ -Zq , the equation (4.12) becomes the homogeneous wave equation as 
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Figure 8. Cross sectional view of the wave front 



(4.16) 


(4.17) 


The solutions of the wave equations in each region must satisfy the boundary conditions, 
which state that the wave and its derivatives must be continuous at all boundaries. Then 
we can find the coefficients of each component wave as they propagate through the 
medium (Figure 8). 




= G'"\ 

dG’’ 

_dG'” 

lz=0 

\z=0 

dz 

z=0 


(4.19) 
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(4.20) 


G'”! 

= g1 , 

dG” 

_aG^ 

\z=L 

\z=L 

az 

7 

Z=L 


We can ignore the coefficient of the left going component at z < -Zq and of the 

right going component at z> L, because the wave is outgoing from the source. In 
addition, by imposing the Sommerfeld radiation condition, the surface integral at infinity 
vanishes [16]. Since the waves are continuous, the Green’s functions and its derivatives 
also have to be continuous, except at the source. Let’s find the derivative of the Green’s 
function. We use equation (4.12) and integrate it on both sides. 

j^ G{v^,v^,z,rQ)dz + \{k^n iz))--{2nv^fG{v^,v^,z,r^,)dz = \Siz +ZQ)dz 

(4.21) 

The left hand side becomes the first derivative of the Green’s function and an integral of 
the Green’s function. The right hand side is the same as the Heaviside Step Function, 
which is a discontinuous function whose value is zero for negative argument and one for 
positive argument. Then the equation becomes 

^G{v^,Vy,z,r^) + \{k^td{z)-{2nvJ^-{2nv^f)G{v^,v^,z,r^)dz = H{z + z^) (4.22) 

The second term in equation (4.22) will vanish, because the Green’s function is 
continuous at any given boundary. But the first term will be continuous as long as the 
boundaries are z < -Zq or z > -Zq . As the first derivative crosses the source where the 

Heaviside Step Function goes from 0 to 1, the difference across will jump to 1. Using the 
same boundaries as the above, we get the following: 

If we let the refractive index obey 

(z<0) 

(0<Z<L) 

Hj (L < z) 

then, we can derive the following equations: 

At z = 0 
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A” +8” =A’^ 


(4.23) 


A\lkX-\F\ -B^4kX-\F\ =A-^k^n^^-F -B^4k^n^^-\F\ (4.24) 


■ 1 1 ■ 

~ A^~ 

■ 1 1 1 " 

-Fk_ 

B'’ 



A" 

B" 


(4.25) 


where = k - F and = k nt -|F| 


m m 


At z = L 


A’^e 


L-ik^nl-\Ff 


+ B"’e 


-iLjk\l-\Ff .TjL^k\F\Ff 


= A‘e 


(4.26) 


im ;,.2„2 I Z7|2 JL^k^nl-\F\^ 


A'”Jk^nt - F e“ 


-B'"Jk^nt - F e 


\T /;,2„2 I 171^ |F| 


= A‘ k^n^r-\F\ e" 


(4.27) 



-iBjFZ 

e ^ 

~A"‘~ 



Rm 


M 3 

~A^~ 


0 


(4.28) 


where Nj. =k^nl - \F^ and 
At z = -Zq 

^iZo<lk^nl-\F\ _ ^-izo'jk^nl-\F\ j^b^iZo^k'‘'nl-\F\ 29) 

We need to use the jump condition, which is the second one in equation (4.18) to obtain 


A\kW-lFre 


2 -izo^Jk^d-lpf rih 


- 2„2 _| p |2 


-B\ k^nt- F e‘ 




2»,2 _\p\^ 


(4.30) 


= -B\ k^nt- F e 


+ 1 
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T 

1_ 

1 - 

'a"' 

T 

_1. 

1 

.1_ 

B’’ 






"5®" 


“O' 

1 

1 


0 J 

+ 

1 


(4.31) 


where • Now we can get the solution for the coefficient of the Green’s 


function. Let’s write our results in the simplified form 

A”' 

5” 



'a'’" 


M, 


= M, 

1 

B’’ 

2 


“a”' 

r 

M 3 

5” 

= M, 


A 

0 




L 

J 

L 

~ 

1 





“O' 


= Mg 


+ 


Bb 


0 


1 


(4.32) 


(4.33) 


(4.34) 


where the reflection and the transmission amplitudes can be expressed in terms of the 
amplitude of the incident wave A* 


'a'’" 



Bb. 


0 

'p 


(4.35) 


b'’=p,,a^ 

(4.36) 


(4.37) 

A,i A,i 


Doing some algebra, we get the reflectance and the reflection coefficient R = r*r . 
where 


5" ( 

Ar+4K)[M-4K) + {-'Pk + ^ 



m j 



\l V 

^)( 

yl^b ■'"V 

'n ] 

l + l 

^ "*"V 


yj ~ Aj 

In ' 

j 
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By the same technique, we can obtain the transmittance and the transmission coefficient 
T = t*t. 


where 


{4K+4K ',)(^+ 


(4.39) 


We also get the other coefficients inside the medium 


B"\ ' ■ '[ 0 


(4.40) 



A™ =P..A^ =^A'' B'”=P2,A^ 

p 

= ^a" 


’ A,i 

A,1 

A” _ 




+A^)(^/A++(-xK+V^)( 

[4y-4yy"''^- 

B” 

2y^4'iy-'iy) 



+x/^)( A+xMi)+(-xK+V^)! 



(4.41) 


(4.42) 


(4.43) 


At the source, we obtain 


^0 =mA[m,m;^m,m,^m,) 


) ~ 1 


(4.44) 


and if we let 




(4.45) 


S=M, 


rol e‘^^yi+4iy) 


^-izo^y 


(i+V^) 


(4.46) 
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then equation (4.44) becomes 


5* 

0 



S 


(4.47) 


where 




A,i 


V ^1.1 y 


\i and 0 = X^^A^-S^^ 


(4.48) 


We know that the refractive indexes outside of the medium are the same: = Nj = N. 

What we need for our simulation are the amplitudes of the incident wave (A*), the 
reflected wave (5*), and the transmitted wave (A^). 




2Viv 


B'’ =r*A'’ 


\ \J m ^ j \ yj m ^ j 


(4.49) 


(4.50) 


A^ =t^A'’ 



(4.51) 


Since the free space Green’s function is a solution of the Helmholtz equation and satisfies 
the radiation condition at infinity, this should be the solution in that region. The following 
are the Green’s functions we are interested in : 


For the region -Zq < z<0. 


G(v^,v^,z,rJ = < 


- -e ^ ' 


24n 


right going 

iZo^N 


(4.52) 


+ 4 n \ -u/i^-ViV e 




left going 


For the region z> L, 
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G{v^,Vy,z,r^) = 




-iL-jN 


N + JN^ e 


24n 


e‘"° ■ 

e'^'‘ ' ' right going 


(4.53) 


We can also get the coefficient A” and B”’ by solving this 


A” 

5” 




5" 


(4.54) 


For the region 0 < z < L, we get these coefficients 


2V^(V^ + 7^) 

2^ + 


5” 



Finally, we obtain the Green’s function in the media 


G(v^,v^,z,ro) = A'”e 




right going wave 


Giv^,v^,z,r^) = B'’'e left going wave 


(4.55) 


(4.56) 


(4.57) 

(4.58) 
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V. THE ANGULAR SPECTRUM 


A. TRANSFORM SPHERICAL WAVES 

The angular spectrum method is a technique for modeling the propagation of a 
wave field. This method involves expanding a complex wave field into a summation of 
many plane waves. As the spherical waves propagate through a wall, they are 
appropriately described by a planer expansion, because waves propagate in all directions 
spherically, and the boundary is plane. We have to transform the spherical waves to plane 
waves, since the reflected wave from a plane media with a spherical wave is also 
spherical [21]. To obtain the angular spectrum g(x,y,z,rg) we need to take the Inverse 

Fourier Transform of the right going wave in the region < z < 0 in equation (4.52). 
We obtain 


g(x,y,z,r^,)=\ f G(v/^,v/ 

J -00 J -oo 


(5.1) 


In this region, the refractive index is n = l, and equation (5.1) becomes 

iZo 


/•OO poo 

g(x,y,z,0 = j_^ J 


izJk^-{I ttv f'-(I ttv f' i 27 r{v^x-\-v y) j j 

? ^ 'e ' dvdv^. 


2Jk^ -{2nvy -ilnvY 


(5.2) 


ik{zQ+z)-j\-{27rv^lkf'-{ln;v^lkf' 

g(x, y, z,ro) = (5.3) 


2kjl-(27rv^/k) -(27rVy/k) 


— 2tc 

Let’s suppose a wave vector k has magnitude and direction cosines (c!;,y0,7), (See 
Figure 9) then k is expressed as 


k =k^x + k y + k z = k{ax + /3y + /z) 


2n 

T 


(ax + /3y + ^Jl-a^-/3^z) 


(5.4) 

(5.5) 
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since direction cosines are «^ + + 7 ^ =1 where 


= k ■ X = ka 
ky=k-y = k/3 
k =k ■z = ky 


Then we can get the following identities 


k^ = ka, k =k/3 ^ a = 


2n:v^ 






Using these identities, equation (5.3) becomes 


/•OO l»00 

g(«,y0,Z,ro) = J ^ J 


^*(Zo+z)Vl-«^-yS^ 

_ ^ik(ax+fiy) 

Ik^X-a^-p^ 

.00 ^if^r(Za+z)^ik(ax+/ 3 y) 


(k ^ 

da 

( k ^ 



ylnj 


\ln J 


d/3 




/•OO fiO 

J -00 J -< 


r 


dadp 


From H. Weyl’s expansion of the spherical waves into planar waves [22] 



Figure 9. The wave vector k 


(5.6) 
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(5.7) 


iK\r-r . 1 

^ f“ r* 1 ^ k[a(x-x')+l3(y-y')+Y(z-z)\) 

|r —r^l j-<k y 

where r = xx + yy + zz, r' = x'x + yy + z'z 
Substituting equation (5.6) into equation (5.7), we get 


dadp 


g(x,y,z,rf^) = -^ 

i47r \r-r^^\ 

and applying the far-field approximation to equation (5.8), we obtain 


g{x,y,z,r^)n 


-ikr-t 


i4nr 


(5.8) 


(5.9) 


B, REFLECTIVITY OF MEDIA 


In the previous chapter, we obtained planer expansions of Green’s functions. To 
calculate the reflectivity on the surface, we need to employ spherical coordinates. In 
equation (4.52), we see the reflectivity {R ) from left going waves as 


R = 




(V^+Viv( 



(5.10) 


where. 


4T^=^ny-{2nv,f-{2nv^f , 4n = -{2nv^f-{2nv^f (5.11) 


The spherical expansion of the reflected wave can be obtained by converting R into 
spherical coordinates. Using the relationship in equation (5.4) and the spherical 
coordinate representation, we obtain 

2nv^=k?,m6cos(/), 2nv^, =ksm6co?,(/), 2nv^=kcos6 (5.12) 


Substituting equation (5.12) into equation (5.11), we get 
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(5.13) 


yjN^ = k^Jnl^ -sin^ 0 , s/n = A:^/l-sin^ 0 = k cos 0 


And inserting these results into equation (5.10), we get the reflectivity as 


R{0) = 


(l-nl) + (nl-l)e 


'2kL\fn 


-sin^ 6* + Vl-sin^ -sin^ ^ - Vl-sin^ e‘ 

/I 2\ .( 2 1 \ i2kLJnlj-sm^ 0 

_ (l-n„) + (n^-lje _ 

i2kLylnf^-sin^ 6 


'2kLyjnf^ -sin^ 9 






nl -sin 0 +cosd\ -\Jnl-sin 0-cos0\ e 


(5.14) 


(5.15) 


C. EVANESCENT WAVES 


The electric field does not abruptly vanish at the boundary, because the phase 
difference between the incident and reflected waves prevents the complete destructive 
interference required to eliminate the transmitted wave. The transmitted field that extends 
beyond the boundary of the dielectric, when a wave is totally internally reflected, is 
known as the evanescent wave [19]. 

If we choose a coordinate system in which the boundary lies in the x- y plane 
and k lies in the x-z plane and suppose we have a transmitted wave of the form 


E, = E.^e 


iikf-r-cot) 


then, we can rewrite this with Snell’s law as 


k,-r= k,x sin 0,+k,z cos 0, 


= k. 


X sin 0. . /sin 6*,. ~ 

- '- + iz. ^^-1 


V 


n 


n 


) 


(5.16) 


(5.17) 

(5.18) 


Substituting equation (5.18) into equation (5.16), we get the transmitted field 

= (5.19) 


40 



where k’ is a constant and k' = k, sin 0,. This wave propagates along the boundary with 
propagation constant k ', diminishing exponentially with constant r. 

Let’s look at the region -Zq < z< L and examine how much the wave evanesces. 
We need equation (4.52) and equation (4.57) for the reflection and transmission. 


g(x, y, z, r„) = 


(5.20) 


g(x, y, z, ro) = 


(5.21) 


For the reflected wave in equation (5.20), there is an angle of incidence that results in a 
transmission angle that is parallel to the surface. If the incident angle increases over the 
critical angle, the transmitted wave will evanesce. In the region z<0, the intensity of 
reflection wave is 


T • 

^=g g 


A Ae 


2(z+Za)^('l!iv^xf-(27iv^yf-k^ 


(5.22) 


As z goes negative, the intensity becomes smaller, so the evanescing wave will be much 
smaller. Therefore in this model system, we will ignore the evanescent waves as 
unimportant. 


h 



Figure 10. The wave vector k^, k^ and k^ 
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Figure 11. Reflection (polystyrene, quartz, glass) 


D. MATERIAL ANALYSIS 


Typical materials of common walls are wood, tile, and concrete. These materials 
have their own refractivity defined by the index of refraction 


c sin 6’. 

n = — = - - 

V sin 0^ 


(5.23) 


The refractive indexes of materials depend on the frequency with which it is measured. 
Unfortunately, indices of refraction in the terahertz region have not been well tabulated. 

Not all light striking a transparent material is refracted as shown in Figure 10 [23]. The 
reflectance, R is related to the index of refraction by 
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R = 


'-k+l 


V^t+1 


zUl 


(5.24) 


(when the incident angle 6. =0). 

From equation (5.24), we can predict that the higher-n materials are the more refractive. 
In this section, we will analyze materials using MATLAB simulations based on equation 
(5.15) and find appropriate frequencies that can penetrate wall materials without much 
loss and with good reflection of objects. For this simulation, we assume that we are using 
plane waves and normal incidence. 

One of the most attractive advantages about terahertz frequencies is their ability to 
penetrate many common materials. But the terahertz pulse can’t penetrate materials with 
high water content, so we assume the materials used in this simulation are pre-dried. 
Table II shows index of refraction for common materials in terahertz frequency [4, 8]. 

Figure 11 shows that as we increase the index of refraction, the reflection 
becomes large. In Figure 12, we used refractive indexes of common wall materials. In the 
case of tile, rock, and sand the average reflection is quite high. But if we look at this 
figure closely, we can find the frequencies that have low reflection from these materials. 
These frequencies are at about 0.94THz, 1.44THz, etc. and can be used for Through-the- 
Wall Imaging. 
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Reflection 


Polystyrene foam 

1.01 

Ochroma pyramidale 

1.08 

Lophira lata (wood) 

1.49 

Skin 

1.4 

Body fluid 

3.9 

Sand 

1.67 

Glass 

1.8 

Quartz 

1.5 

Tile 

1.6 

Stone 

1.6 


Table 2. Index of refraction for common materials. 
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Figure 12. Reflection (lophira, tile, sand) 


44 
















We also used the index of refraction of human skin and body fluid. Figure 13 
illustrates the reflection from body fluid is very high at frequency 0.94THz. Consequently, 
if we choose 0.94THz frequency for Through-the-Wall Imaging, we will get more data 
from the human body and get a good image. 
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Figure 13. Reflection (body fluid, skin) 
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VI. INTEGRAL EQUATION 


A. WAVE PROPAGATION EQUATION 

In this thesis we are assuming that the complex dielectric constant and magnetic 
permeability vary slowly over a wavelength of the electromagnetic wave, and we are 
neglecting polarization effects. Then we can write the Helmholtz equation as 

V^E{r) + enl.{r)E{r)=Q ( 6 . 1 ) 

If we add k^E{r) on both sides we get 

V^E{r) + eE{r) = -k\n„,,{rf -l)E{r) (6.2) 

0(r) 

where 0(r) is the target’s object function and is denoted by 

0{r) = k\nl.-l). (6.3) 

Then the total electric field ^(r) solution to equation (6.2) is 

E(.r) = E,„,(r) + g(r -rJO(rJE(rJdr„ (6.4) 

where g is the Green’s function of the Helmholtz equation and is the solution of the 
differential equation 

(W^^)g(r-r„) = -^(r-r„) (6.5) 

where (5’(r-rj) is the Dirac delta function. In our model system, the Green’s function 

varies as the observation point changes. From the previous chapter, we found Green’s 
functions in the spatial frequency domain as Eq.(4.52), Eq.(4.53), Eq.(4.57), and 
Eq.(4.58). So in order to get g(r-rQ), we have to transform these Green’s functions 
using the Inverse Eourier Transform as 

g(X-^o) = \\\G(y,,Vy,z,r^)e^’‘^'''"^''^"^dv^dVy ( 6 . 6 ) 
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We denote each Green’s function by their domain as 

^/(r-ro) for z<0 
^//(r-r#) for 0<z<L 
^///(r-ro) for L<z 

and we will use gjj, (r - r,, ) for our model system in the following section. 

B, THE FIRST BORN APPROXIMATION 

The Born approximation is perhaps the simplest and most widely used EM 
scattering approximation. It is based on the assumption that scattered electric fields inside 
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scatterers are negligible compared to the normal/background electric fields. Hence the 
second term on the right hand side of equation (6.4) is small in comparison to the first 
term within the integral. In the integrand can approximate 

E{r)^E,Jr) (6.7) 

Therefore the total electric field in equation (6.4) is calculated as 

E(y) = (r) + g(r-r^)0(r^)E.^^(r^)d\ (6.8) 

Consequently, the scattered field, E ^, can be identified with 

=£^(r-r„)6>(r„)£.„^(r„)dV„ (6.9) 

C. THE FIRST RYTOV APPROXIMATION 


The Rytov approximation assumes that the incident wave perturbation caused by 
the target can be described by a change of phase in the reference wave. The Rytov 
approximation is obtained by representing the total field as a complex phase [22] 


^(r) = e 




( 6 . 10 ) 


We consider the scattering of a scalar wave in an inhomogeneous medium. The field 
satisfies the reduced wave equation 


V"E(r) + k"E(r) = 0 

Substituting equation (6.10) into equation (6.11), we get 


Using the identity 


vV^''^+kV“'^ =0 




( 6 . 11 ) 


( 6 . 12 ) 


(6.13) 


we substitute equation (6.13) into equation (6.12) and divide by on both sides and 
get the following Riccati equation 


vV + (V^)(V^) + k" =0 


(6.14) 


We denote the refractive index as 


k{r) = kijiir) = k^il + ndir)) 


(6.15) 
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The zero perturbation field for n(r) = 1 can be written as E.(r) = , and its phase will 

satisfy equation (6.14) 

V^^.+(V(^f+k^=0 (6.16) 

We can express the total complex phase, (j ), as the sum of the incident phase and the 
scattered complex phase (j)^ 

(/> = (/>,+(!>, (6.17) 

Now substitute equation (6.17) into equation (6.14) 

+ (6.18) 

Expanding this equation 

vV,+vV,+(V^zi,)"+(V^zijV2(V^zi,)(V^ziJ + ^^=0 (6.19) 

and subtract equation (6.16) from equation (6.19), we get 

VV, +2(V(zi,)(V(ziJ + (^2 -^„^) = -(V(ziJ^ (6.20) 

VV, + 2(V(zi,)(V(ziJ = -1)) (6.21) 

Using the identity 

V2(E,^ziJ = + 2(VE,)(V^ziJ + E,VV, (6.22) 

and assuming a plane wave for the incident field, so that 

E. = (6.23) 

then we find the following expression 

V^E.=-kf,E. and VE.=E.V(/>. (6.24) 

So equation (6.22) becomes 

2E,V<P,-V<P, + E,V^<P^=V\E,<P^) + e,E,<P^ (6.25) 

Now substituting equation (6.25) into equation (6.21), we get 

V\E,<PJ + k\E,<PJ = -E, {(V<Pf + k^,(n^ -1)} (6.26) 

Solving this equation using the free space Green’s function, we get 

£,(r¥,(r) = j«(r-r.)£,(r,){(V(«,(r.)"+*:„"(«=(r.)-l))*. (6.27) 

£'.(r)-' £',.(r)-' 
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By the Rytov approximation, the first term in equation (6.28) is very small, so equation 
(6.28) can be approximated as 

«<'■ - >■.)£,-D*. (6.29) 

If we denote the target object function 0(r) as 

0(r) = e^(n\rj-l) (6.30) 

Then equation (6.29) becomes 

f ^(r - r# )E, (r„)6>(r„)Jr„ (6.31) 

The corresponding solution for the total field is 

E(r) = E.(r)e^- (6.32) 

But it is only valid when (V (j )^)^ □ kl (n^(r^) -1) 

From equation (6.10), the total electric field is given 

E{y) = E.{y) + E^r) = (6.33) 


Rearranging this equation for the scattered field, E^{r ), we get 


£,(r) = £,(r)(^’^'«-l) 


(6.34) 


Substituting equation (6.31) into equation (6.34), finally we get the scattered field as 


EM) = E,{v,) 


e 


1 


1 


«(r-ro)£,(io)0(ro)rfr„ 



(6.35) 


Note that when the argument of the exponential is small, then we can approximate 

r'^j«(r-ro)£,(>'o)0(ro>(r„ ^ 




l + TrTT f ^(r-ro)£',(r„)6>(r„)Jr„ 
EXr)^ 


(6.36) 


and, in this case, the Rytov approximation reduces to the first Born approximation. 
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VII. SIMULATION 


A. RECONSTRUCTION ALGORITHM 

For many reasons, it is very useful to convert functions in the time domain into 
functions in the frequency domain. The primary reason is that it gives us an easy way to 
solve difficult problems. For instance, when we meet the convolution integral in the time 
domain, we can express this as a simple multiplication integral in the frequency domain. 

To apply the Rytov approximation for object reconstruction, we need to estimate 
the scattered field , which is denoted in equation (6.35). From equation (6.35) we 
have 

£,(r)ln 1^+1 =f *(r-r,)£,(r,)0(r,>Jr, (7.1) 

LE,.(ro) J 

Since a kernel is formed as a convolution between the Green’s function, the scattering 
potential, and the incident electric field, we make a reconstruction algorithm in the 
frequency domain by Fourier Transform of the measurement data. The left term in 
equation (7.1) is our measurement data. If we denote this by m(r), and then 

m(r) = f g(r -rJE.(rJO(rJdr„ (7.2) 

We extend the volume integral to infinity and take 2D Fourier Transform by assuming 
receivers will be place at fixed points Then equation (7.2) becomes 

M(v„k^) = G(v^,i/^)X(v^,v^) (7.3) 

where, 

00 

= j J E,(x„ y,)0(x„ y,)e-‘^^^''^^^''^^^dxdy (7.4) 

—00 

By using equation (7.3) and applying equation (2.9) we can get X(v^,v^) as 
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) = |(G(v^, v^, fG(v^,v^))' G(v^v^,) 


(7.5) 


Substituting equation (7.5) into equation (7.4), we can get the “best object” 0(rg) and 
taking the Inverse Fourier Transform of this as 


OW = j [ l(G(v,.y,fG(v,.vjy‘ G(v,.vj}M(v,.yJ 


i2;T(y^x+y y) 


dxdy (7.6) 


Finally, we get the “best object” Oir^) in the spatial domain. The incident wave 
£',(r) = Aebecomes E[{r) = after passing through a transparent media with 

transmittance t. The Green’s function to be used in equation (7.6) is 


G(v^,v ,z,ro) = 




il4n 


N+jN. 




-dN 


24n 




(7.7) 


B. MATLAB CODE IMPLEMENTATION 

The total field at the receivers will be composed of the incident field from 
transmitters, the reflected field from the wall, the scattered field from the target, and the 
evanescing field. In the previous chapter, we showed the evanescent field is negligible in 
comparison with the incident field. 

For our image processing, we are only interested in the scattered field data and 
use equation (6.35) and equation (7.6) to create scattered field images. Due to the 
limitation of computer ability, we are restricted in the number of scattering points to the 
sampling frequency, and bandwidth to be examined. 
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(b) 

Figure 15. Amplitude of field scattered from 5 points using Born (a) and Rytov (b) 

approximations 

Figure 15 shows 5 scattering points applying the Born and Rytov approximations 
at z = 2m . We will add Gaussian noise to this data and use Tikhonov regularization, 
truncated SVD methods, and the L-curve method discussed in chapter 2. All MATLAB 
codes are presented in Appendix D. 
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c. 


SIMULATION RESULTS 


To compare the Born and Rytov approximations, first we simulate with equation 
(6.9) and equation (6.35) separately. We get Figure 15 (a) from the Born approximation 
and Figure 15 (b) from the Rytov approximation. These two figures look very similar, 
and they localize the scatterers equally well. But the strength of each peak differs by 
1.2% in both figures. The difference is formed by subtracting equation (6.35) from 
equation (6.9) and simulating again. Figure 16 shows us the differences between them. 
The differences of each peak are approximately 0.16%, and both approximations perform 
equally well. 

D. SIGNAL TO NOISE RATIO 

For the simulation of measuring a scattered electric field, we have to add noise to 
the original signal from the Rytov approximation, since signals are corrupted by noise in 
nature. We will use the Gaussian noise that is commonly used as an additive noise model 
in physics. Since the Signal to Noise Ratio (SNR) is defined as the ratio of signal power 
to noise power, we calculate the SNR by taking the ratio of the sum of the signal and the 
Gaussian noise in the frequency domain. Furthermore we find the minimum SNR to 
recognize the image of objects by increasing a noise factor. 
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Figure 16. Differences between Born and Rytov approximations 


Noisy Data with SNR= -20.0 dB 



Figure 17. Scattering points with Gaussian noise SNR=-20dB 
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E. 


TIKHONOV REGULARIZATION 


Figure 17 shows the scattered signals with noise. We use Tikhonov regularization 
to mitigate these noise effects. During this process there is a trade-off between the 
regularized output and the original data. In order to optimize the trade-off, we have to 
choose the best regularization parameters. To find the regularization parameters we will 
use L-curve method. The L-curve method is based on a log-log plot of the norm of the 
residual versus the solution for a range of values of regularization parameters [14]. The 
common expression of the Tikhonov regularization is 

F^=\\AX-b\l+a\\x\l (7.8) 


II l|2 II ||2 

where ||AA-£>||2 is a residual norm, is a regularized solution, and a is the 

regularization parameter. The “horizontal” part of the L-curve is dominated by 
regularization errors that occur from over-smoothing and the “vertical” part is dominated 
by perturbation errors that occur from under-smoothing. Hence, we choose our optimal 
solution as the one at the corner of the L-curve in Figure 18(b). So a particular 
regularization parameter is a = 4.6. Using this parameter we can get a good recovered 
image in Figure 19. However, Figure 20 shows us that if we use regularization parameter 
smaller than the optimal one, we will get an unstable image. 
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Figure 18. 


(b) Residual norm 

Regularization parameter (a) and L-curve (b) 
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F. TRUNCATED SVD 

Since the operator A of TWI has large singular values, we can apply the 
truncated SVD method. Using the singular value spectrum plot in Figure 21, we can find 
values /C = 0.34 and ,/C = 0.10 where the singular values change rapidly, and truncate 
with these values. Figure 22 is obtained by truncating /C = 0.34 ,and Figure 23 is 
obtained by truncating /C = 0.10 . These two figures are good examples of optimal 
truncation parameters. 

On the contrary. Figure 24, obtained at /C = 100 , shows that a too-large /C affects 
the accuracy of estimation. 
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Figure 21. Singular value spectrum of transfer operator 
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TSVD w/ truncation at l<appa= 0.34 
SNR=-20.0 dB 
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Figure 22. TSVD K =0.34 at SNR=-20dB 


TSVD w/ truncation at kappa= 0.10 
SNR=-20.0 dB 
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Figure 23. TSVD K =0.10 at SNR=-20dB 




TSVD w/ truncation at kappa= 100.00 
SNR=-20.0 dB 
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Figure 24. TSVD K =100 at SNR=-20dB 
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VIII. CONCLUSION 


In recent years, a variety of techniques for THz imaging have been developed for 
the detection of weapons, explosives, and drugs at security check points and for the 
detection of illegal immigrants at the border. This thesis is a review of some of this work. 
We think our model may be useful for military and security officers to detect enemies 
hiding in a room or behind a wall. 

In this thesis we have applied the Bom and Rytov approximations to the problem 
of imaging from scattered field measurements. Also we examined the Tikhonov 
regularization parameter selection with the L-curve method. The simulation result 
obtained by the Rytov approximation has been compared to one yielded by the Bom 
approximation. The difference between them is insignificant. Therefore the Rytov 
approximation appears to be no better than the Born approximation. It would be useful to 
apply this simulation to real data in the future. In general, the Rytov approximation is 
more accurate than the Bom approximation especially, at higher frequencies [24]. 

Further modeling efforts are also required. We have assumed the medium (wall) 
is linear, nonmagnetic, and nonconductive. But in reality, wall materials are not exactly 
like this. Some materials are nonlinear, magnetic, and conductive. So for this case, the 
wave equation has to be adjusted to derive a new Green’s function to fit this model, and 
we have to apply vector theory instead of scalar theory to the wave equation. 

For more accurate image reconstmction, we also need to consider dispersive 
effects of different kinds of walls and have to include multiple scattering effects to the 
Lippman-Schwinger equation. 
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APPENDIX A. SINGULAR VALUE DECOMPOSITION 


Let K e □ is a linear operator which maps vectors / e □ ^ to vectors d & . 

We will consider two symmetric matrices ¥J¥^{NxN) and KK^(MxM). Since these 
matrices are square and symmetric, we can get eigenvectors and eigenvalues of each. If 
/I, is the eigenvalue and f. is the eigenvector of K^K . 

The eigen equation for K^K is 

(A.l) 

Applying the same method, we get the eigen equation for KK^ 

¥^¥Jd. = y.d. (A.2) 

If we take the eigenvector of K^K, and the eigenvalue A- then we get vector 
Kf. 7^0 , this can be proved here 

(KK")K/;. =K(K"K)/;. =Ka,)/;. =i,(K/;.) (A.3) 

Hence Kf. is an eigenvector of KK^ and each eingenvalue A. of K^K must be an 
eigenvalue of KK^ . Similarly, this is also true each eigenvalue y. of KK^ must be 
an eigenvalue of K^K . 

If there are r non-zero eigenvalues of K^K (or KK^), where r<m and r<n, we can 
express it like this 

^1 =/i 

A+i = 0 y,^i = 0 

^„=o r„=o 
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Therefore the linear equation is expressed 


/ = YJd d = Kf 


(A.4) 


Then we normalize eigenvectors 


K^d 

K^d 


d = 


K/ 


The singular values cr^ of K are defined by 

||ka I = \\ = (^k = ^^\ =yjn 

So equation (A.5) becomes for K < M, A . 

K/, =o-,d, 

but for K > M, A , the equation becomes 

K/^ = 0 for k = r + l,---N 
'KJd^=0 for k = r+ 


(A.5) 

(A.6) 

(A.7) 

(A.8) 

(A.9) 


On equation (A.7), we take K/^ =(^k^k multiply by // on both sides. Then it 
becomes 


K"=Zct,/,< 

/=! 

And the other one in equation (A.7) becomes 

i=\ 

We can write equation (A. 11) in matrix form as 





0 

o' 

-... ...- 

K = 


0 

^2 

0 

... ... 



0 

0 




(A. 10) 


(A. 11) 


(A. 12) 


or K = DUF^ 

where, t/ is a diagonal matrix of the singular values, and D and F consist of the 
eigenvectors corresponding to those singular values. / and d are the orthonomal vectors 
and called right singular vectors and left singular vectors respectively. 
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We can make the inverse of K as 


K ‘ = (D[/F0 ‘ 


K‘=(D[/F0‘ (A.13) 

We know and D are orthogonal matrices, because they are composed of 
eigenvectors which are orthogonal. Equation (A.13) becomes 

K ' =(E^) 't/ 'D ' (A.14) 

K' = Ft/'D^ (A. 15) 

is a diagonal matrix composed of —. Since singular values of K are in the order of 


cr 


CT; >0-2 >--->cr^ >0 


(A. 16) 


elements in a matrix U ' will increase infinitely [10]. Hence, the norm of inverse 
mapping also will increase. This is why we need a regularization method for the 
inverse imaging. 
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APPENDIX B. GREEN’S FUNCTION 


Green’s function satisfies a wave equation driven by a point source. 


+k^)i//ir) = Sir) 

(B.l) 

To obtain a solution to wave equation need to seek the Green’s 

solution to the following equation. 

function that is the 

(V^+k^)G(r,r') = -d(r-r') 

(B.2) 

If we let arbitrary source Sir) 


Sir) = ^ Sir')dir -r')dr' 

(B.3) 

Then we can write 


i//ir) = -J Gir,r')Sir')dr' 

(B.4) 

To find a solution of equation (B.2), solve it in spherical coordinates with the origin at r' 


iV^+k^)Gir) = -dir) = -dix)diy)diz) (B.5) 


Due to the spherical symmetry of a point source, G(r) is also spherical symmetric. 

In spherical coordinates equation (B.2) becomes as 

-^irG) + eG = -dir) (B.6) 

r dr 

Everywhere, except at r = 0, the equation is the same as the homogeneous equation 

^irG) + k\rG) = Q (B.7) 

dr 

So the solution is shown below 

ikr 

Gir)=A — + B - (B.8) 

r r 

Since sources are absent at infinity, only an outgoing solution can exist. 

ikr 

Gir) = A— (B.9) 

r 

To find the constant A , substitute equation (B.9) into equation (B.2) and integrate 
equation (B.2) over a small volume. 
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r 


(B.IO) 


Since dV = Anr^dr , the second integral vanishes and the first integral is converted into a 
surface integral using Gauss’ Theorem 

2 d Ae'"' 1 

hm4;rr^-= -l (B.ll) 

'•^0 dr r 


Finally we get A = 

Therefore, the Greens’ function G is 

ik\r-r'\ 

G(r-r')=^ -H (B-12) 

4;r|r-r I 

Consequently, the final solution to equation (B.l) is 

ik\r-r'\ 

¥{r) = - f - yS{r')dr' (B. 13) 

* 4;r r - r 
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APPENDIX C. FOURIER TRANSFORM IN 2-D 


When dealing with both linear and nonlinear systems, it is useful to break down a 
complicated input into more simple inputs. This process is well known as the Fourier 
Transform. The Fourier Transform of a complex-valued function g of two independent 
variables x and y is defined by 

00 

•^{^} = J J six, 


Similarly, the Inverse Fourier Transform of a is defined as 

00 

{G} = J j (C.2) 

-oo 

But all functions cannot be transformed from the spatial domain into Fourier domain. In 
order to be transformed, the function should satisfy “existence conditions” [20]. 

1. ‘ g ’ must be absolutely integrable over the infinite (x, y) plane. 

2. ‘ g ’ must have a finite number of discontinuities and a finite number of maxima 
and minima in any finite rectangle. 

3. ‘ g ’ must have no infinite discontinuities. 

There are a few useful properties of the Fourier Transform. These will provide the basic 
idea for the manipulation of the Fourier Transforms and can save time solving Fourier 
problems. 

1. Linear Theorem : 

T{ag+ldk) = aT{g)^ldT{k) (C.3) 


2. Similarity Theorem 

If J^{g(x,y)} = G(i/^,i/ ),then T{giax,by)] = ^G(^Xf-) 

\ab\ a b 


(C.4) 
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3. 


Shift Theorem : 


If T{g{x,y)] = G(y^,Vy),t^&ri T{g{x-a,y-b)] = G(y^,Vy)e 


-i27r{v^a+v b) 


(C.5) 


Parseval’s Theorem 


CJU 

If {g(x,y)} = G(v^,v^) j J \g(x,y)fdxdy = j ^ \G(v^,v^f\dvJv^ (C.6) 


5. Convolution Theorem : 

If J^{g{x,y)] = G{v^,v^) and j)} = , then 



(C.7) 


6. Autocorrelation Theorem 

If J^{g{x, j)} = G(i/^,i/^),then 


j j g{a,b)g*{a-x,b-y)dadb\ = \G(,v^,v^) 


similarly, y)|^| = J ^\G{a,bf'\^G*{a-v^,b-Vy)dadb 

-00 


(C.8) 

(C.9) 
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APPENDIX D. MATLAB CODES 


%%%%%%%%%%%% Matlab code for material reflection %%%%%%%%%%%%% 
clear; 

% Calculating the Reflection for plane waves in various surface 
% Equations used from Chapter 5 

% Creating various index of refraction for different materials 
% Polystyrene foam, Quartz, Glass 
nm=[1.01 1.5 1.8]; 

% Lophira alata(wood). Tile, Sand 
%nm=[1.49 1.6 1.67]; 

% Skin, Body fluid 
% nm=[1.4 3.9]; 

% Thickness of the material 
L=0.5; 

c=3.0e8; % Speed of light 

figure(l) 

nu=linspace( lell,5el2,l 00); 
k=2*pi*nu./c; 

num=( 1 -nm( 1 )^2)+(nm( 1 )^2-1) *exp(i*2*k*L*nm( 1)); 
den=(nm( 1)+1 )^2-((nm( 1)-1 )^2)*exp(i*2*k*L*nm( 1)); 
r=num./den; 

R=r.*conj(r); 
hold on 
plot(nu,R,'-'); 

gtext(sprintf('n=% 2.2f ,nm( 1))); 

num=(l-nm(2)^2)+(nm(2)^2-l)*exp(i*2*k*L*nm(2)); 
den=(nm(2)+1 )^2-((nm(2)-1 )^2)*exp(i*2*k*L*nm(2)); 
r=num./den; 

R=r.*conj(r); 
hold on 
plot(nu,R,'-.'); 

gtext(sprintf('n=%2.2f,nm(2))); 

num=( 1 -nm(3)^2)+(nm(3)^2-1 )*exp(i*2*k*L*nm(3)); 
den=(nm(3)+1 )^2-((nm(3)-1 )^2)*exp(i*2*k*L*nm(3)); 
r=num./den; 

R=r.*conj(r); 
hold on 
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plot(nu,R,':'); 

gtext(sprintf('n=%2.2f,nm(3))); 
hold off 

xlabel('frequency (\nu) Hz'); 
ylabel('Reflection'); 

title('Frequency vs Reflection from the dielectric') 
gtext(sprintf('L = %1.2f,L)); 

%%%%%%%%%%% Matlab code for Bom and Rytov approximations %%%%%%%% 


clear; 

%parameters 

nu=940e9; 

lambda=3.0e8/nu; 

%k=2*pi/lambda; 

k=3000; 

n=1.67; 

L=0.5; 

Z=2; 

A=l; 

signal=0; 

noise=0; 


% Frequency of the wave propagation. 

% In terms of wave length. 

% Magnitude of the propagation vector. 

% Material index of refraction. 

% Thickness of the wall. 

% Measurement plane at Z=constant. Adjust accordingly 
% depending on position of the receiver. 

% Amplitude of the incident plane wave. 


% Defining Coordinate Systems in the spatial frequency domain. 
ii=10; % Spatial frequencies for meshgrid input. 

fb=ii; % Frequency band. 

%fn=fb; 

fn=2*fb; % Nyquist rate. 

j=l/fn; % Separation must be at least this distance to prevent aliasing. 

[fx,fy]=meshgrid(-ii:j:ii,-ii:j:ii); % Creating coordinate system. 

% Defining the Green's function. 

x0=0; y0=0; z0=-3; % Set this position accordingly to the problem. 

N=sqrt(k^2-(2*pi*fx).^2-(2*pi*fy).^2); % Look in chapter 4 for reference. 

Nm=sqrt((n*k)^2-(2*pi*fx).^2-(2*pi*fy).^2); % Look in chapter 4 for reference. 
num=4*N.*Nm*exp(-i*N*L); % Numerator of the equation. 

den=(Nm+N).^2.*exp(-i*Nm*L)-(Nm-N).^2.*exp(i*Nm*L); % Denominator of the 

equation. 

% Required terms for taking Fourier transform 
A0=exp(-i*2*pi*fx*x0).*exp(-i*2*pi*fy*y0); 

G=(num./den).*exp(i*N*Z).*exp(i*N*z0)./(2*N).*A0; % Transmission Green's function. 
%% Synthetic Data Generation 
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% Transmittance for a plane wave through a given media 
t=4*n*exp(-i*k*L)./((n+1 )^2*exp(-i*n*k*L)-(n-1 )'^2*exp(i*n*k*L)); 

% scattering Points in the spatial frequency domain. 

xl=-l; yl=-l; zl=-l; 

x2=-1; y2= 1; z2=- 3; 

x3=0; y3=0; z3=-5; 

x4=l; y4=l; z4=-1.3; 

x5=l; y5=-l; z5=-2.2; 

% Incident field after the penetration of the media. 

El=t*A*exp(-i*k*zl); 

E2=t* A*exp(-i*k*z2); 

E3=t* A*exp(-i*k*z3); 

E4=t* A*exp(-i*k*z4); 

E5=t*A*exp(-i*k*z5); 

% Position of the scatter points in frequency domain for the given 
% measurement plane at Z=constant 

% Since the scatter points are in terms of dirac delta functions, 

% the field measurement would become the summation of the Green's function 
% and the incident wave after performing the Rytov Approximation 
G1 =(num./den). *exp(i*N*Z). *exp(i*N*z 1) ./(2*N); 
G2=(num./den).*exp(i*N*Z).*exp(i*N*z2)./(2*N); 
G3=(num./den).*exp(i*N*Z).*exp(i*N*z3)./(2*N); 

G4=(num./den). *exp(i*N*Z). *exp(i*N*z4) ./(2*N); 
G5=(num./den).*exp(i*N*Z).*exp(i*N*z5)./(2*N); 

% Required terms for taking Eourier transform 

A1 =exp(-i*2*pi*fx*x 1). *exp(-i*2*pi*fy *y 1); 

A2=exp(-i*2*pi*fx*x2).*exp(-i*2*pi*fy*y2); 

A3=exp(-i*2*pi*fx*x3).*exp(-i*2*pi*fy*y3); 

A4=exp(-i*2*pi*fx*x4).*exp(-i*2*pi*fy*y4); 

A5=exp(-i*2*pi*fx*x5).*exp(-i*2*pi*fy*y5); 

%% Generating Data with Born approximation 

D=(E1.*G1.*A1+E2.*G2.*A2+E3.*G3.*A3+E4.*G4.*A4+E5.*G5.*A5); 

%% Generating Data with Rytov approximation 

DD=E 1 .*(exp(G 1 .* A1)-1 )+E2. *(exp(G2. * A2)-1 )+E3. *(exp(G3. * A3)- 

1 )+E4.*(exp(G4.* A4)-1 )+E5. *(exp(G5. * A5)-1); 

%% Generating differences by subtracting Born with Rytov 
% Generating Data with Bom 
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B=(E1.*G1+E2*G2+E3.*G3+E4.*G4+E5 *G5); 


% Generating Data with Rytov 

R=EE*(exp(Gl)-l)+E2.*(exp(G2)-l)+E3 *(exp(G3)-l)+E4*(exp(G4)- 
l)+E5 *(exp(G5)-l); 

Ext=(B-R).*(Al+A2+A3+A4+A5); % Data with extra terms 

% Setting matrix dimension 
[MM,NN]=size(D); 

D=D(1:MM,1:NN); % Change matrix dimension if warranted. 

% Converting bin numbers to real axis numbers in spatial domain(meters). 

% Range has to be up to Nyquist frequency. 
ax=linspace(-fn/2,fn/2,length(D)); 

%% End of Data Generation 

figure(l) 

dl=ifft2(D); % Data in the spatial domain. 

dl=ifftshift(dl); 

mesh(ax,ax,abs(dl)); % Amplitude of field scattered from 5 points using Bom 
view(-70,20) 

figure(2) 

d2=ifft2(DD); % Data in the spatial domain. 

d2=ifftshift(d2); 

mesh(ax,ax,abs(d2)); % Amplitude of field scattered from 5 points using Rtov 
view(-70,20) 

figure(3) 

d3=ifft2(Ext); % Data in the spatial domain. 
d3=ifftshift(d3); 

mesh(ax,ax,abs(d3)); % Differences between Born and Rytov approximations 
view(-70,20) 

%%%%%%%%%%% Matlab code for regularization mathods %%%%%%%%%%%% 
clear; 

%parameters 

nu=940e9; % Erequency of the wave propagation. 

Iambda=3.0e8/nu; % In terms of wave length. 

%k=2*pi/lambda; % Magnitude of the propagation vector. 
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k=3000; 

n=1.67; % Material index of refraction. 

L=0.5; % Thickness of the wall. 

Z=2; % Measurement plane at Z=constant. Adjust accordingly 

% depending on position of the receiver. 

A=1; % Amplitude of the incident plane wave. 

signal=0; 

noise=0; 


% Defining Coordinate Systems in the spatial frequency domain. 
ii=10; % Spatial frequencies for meshgrid input. 

fb=ii; % Frequency band. 

fn=2*fb; % Nyquist rate. 

j=l/fn; % Separation must be at least this distance to prevent aliasing. 

[fx,fy]=meshgrid(-ii:j:ii,-ii:j:ii); % Creating coordinate system. 

% Defining the Green's function. 

x0=0; y0=0; z0=-3; % Set this position accordingly to the problem. 

N=sqrt(k^2-(2*pi*fx).^2-(2*pi*fy).^2); % Look in chapter 4 for reference. 

Nm=sqrt((n*k)^2-(2*pi*fx).^2-(2*pi*fy).^2); % Look in chapter 4 for reference. 

num=4*N.*Nm*exp(-i*N*L); % Numerator of the equation. 

den=(Nm+N).^2.*exp(-i*Nm*L)-(Nm-N).^2.*exp(i*Nm*L); % Denominator of the 

equation. 

% Required terms for taking Fourier transform 
A0=exp(-i*2*pi*fx*x0).*exp(-i*2*pi*fy*y0); 

G=(num./den).*exp(i*N*Z).*exp(i*N*z0)./(2*N).*A0; % Transmission Green's function. 
% Synthetic Data Generation 

% Transmittance for a plane wave through a given media 
t=4*n*exp(-i*k*L)./((n+1 )^2*exp(-i*n*k*L)-(n-1 )^2*exp(i*n*k*L)); 


% scattering Points in the spatial frequency domain. 

xl=-l; yl=-l; zl=-l; 

x2=-l; y2=l; z2=-3; 

x3=0; y3=0; z3=-5; 

x4=l; y4=l; z4=-1.3; 

x5=l; y5=-l; z5=-2.2; 

% Incident field after the penetration of the media. 
El=t*A*exp(-i*k*zl); 

E2=t* A*exp(-i*k*z2); 

E3=t* A*exp(-i*k*z3); 

E4=t* A*exp(-i*k*z4); 

E5=t*A*exp(-i*k*z5); 
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% Position of the scatter points in frequency domain for the given 
% measurement plane at Z=constant 

% Since the scatter points are in terms of dirac delta functions, 

% the field measurement would become the summation of the Green's function 

% and the incident wave after performing the Rytov Approximation 

G l=(num./den). *exp(i*N*Z). *exp(i*N*z 1 )./(2*N); 

G2=(num./den).*exp(i*N*Z).*exp(i*N*z2)./(2*N); 

G3=(num./den).*exp(i*N*Z).*exp(i*N*z3)./(2*N); 

G4=(num./den).*exp(i*N*Z).*exp(i*N*z4)./(2*N); 

G5=(num./den).*exp(i*N*Z).*exp(i*N*z5)./(2*N); 

% Required terms for taking Fourier transform 

Al=exp(-i*2*pi*fx*xl).*exp(-i*2*pi*fy*yl); 

A2=exp(-i*2*pi*fx*x2).*exp(-i*2*pi*fy*y2); 

A3=exp(-i*2*pi*fx*x3).*exp(-i*2*pi*fy*y3); 

A4=exp(-i*2*pi*fx*x4).*exp(-i*2*pi*fy*y4); 

A5=exp(-i*2*pi*fx*x5).*exp(-i*2*pi*fy*y5); 

% Generating Data 

D=E1. *(exp(G 1. * A1)-1 )+E2 *(exp(G2 * A2)-1 )+E3. *(exp(G3. * A3)- 
1 )+E4 *(exp(G4 *A4)-1 )+E5.*(exp(G5.*A5)-1); 

% Setting matrix dimension 
[MM,NN]=size(D); 

D=D(1:MM,1:NN); % Change matrix dimension if warranted. 

% Converting bin numbers to real axis numbers in spatial domain(meters). 

% Range has to be up to Nyquist frequency. 
ax=linspace(-fn/2,fn/2,length(D)); 

n=(randn(MM,NN))+i*(randn(MM,NN)); % noise in complex values in the frequency 
domain 

% set the signal to noise in dB 
SNR=-20; 

S=abs(D.*conj(D)); 

NO=abs(n.*conj(n)); 

signal=sum(sum(S)); 

noise=sum(sum(NO)); 

factor=sqrt((signal/noise)*10^(-SNR/10)); 

n=n*factor; 

d=D+n; % data with noise in the frequency domain 

% End of (noisy) Data Generation 

figure(2) % Data corrupted with noise in spatial domain. 
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dl=ifft2(d); 

dl=ifftshift(dl); 


% Data in the spatial domain. 


imagesc(ax,ax,abs(dl)); 
colormap (1-gray); 

%grid on 
axis([-2 2-2 2]) 

title(sprintf('Noisy Data with SNR= %3.1f dB',SNR)) 

% Tikhonov Regularization 
[u,s,v]=svd(G); 

alpha=4.6; % Regularization parameter 

Reg=inv(alpha*eye(size(G'*G))-i-G'*G)*G'; 

figure(3) 

dd=Reg.*d; 

ddl=ifft2(dd); % Data in the spatial domain 

ddl=ifftshift(ddl); 

imagesc(ax,ax,abs(ddl)); 

colormap (1-gray); 

%grid on 
axis([-2 2-2 2]) 

title(sprintf('Tikhonov Regularization alpha= %4.1f and SNR=%3.1f dB',alpha,SNR)) 
% Truncated SVD Method 

kappa=0.34; % Parameter for truncation of singular values 

[u,s,v]=svd(G); 

% inverse of the singular value matrix 
sinvlu=zeros(NN,MM); 
for j=l:rank(s) 
if l/s(j,j) < kappa 
sinvlu(j,j)=l/s(j,j); 
else 

sinvlu(j,j)=0; 

end 

end 


td=(v*sinvlu*u').*d; 

TD=ifft2(td); % Data in the spatial domain 

TD=ifftshift(TD); 


figure(4) 

imagesc(ax,ax,abs(TD)); 
colormap (1-gray); 
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axis([-2 2-2 2]) 

title(sprintf('TSVD w/ truncation at kappa= %3.2f \n SNR=%3.1f dB',kappa,SNR)) 

%%%%%%%%%%%%%% Matlab code for L-curve plot %%%%%%%%%%%%%% 


clear; 

%parameters 

nu=940e9; 

lambda=3.0e8/nu; 

%k=2*pi/lambda; 

k=3000; 

n=1.67; 

L=0.5; 

Z=2; 

A=l; 

signal=0; 

noise=0; 


% Frequency of the wave propagation. 

% In terms of wave length. 

% Magnitude of the propagation vector. 

% Material index of refraction. 

% Thickness of the wall. 

% Measurement plane at Z=constant. Adjust accordingly 
% depending on position of the receiver. 

% Amplitude of the incident plane wave. 


% Defining Coordinate Systems in the spatial frequency domain. 
ii=10; % Spatial frequencies for meshgrid input. 

fb=ii; % Frequency band. 

fn=2*fb; % Nyquist rate. 

j=l/fn; % Separation must be at least this distance to prevent aliasing. 

[fx,fy]=meshgrid(-ii:j:ii,-ii:j:ii); % Creating coordinate system. 

% Defining the Green's function. 

x0=0; y0=0; z0=-3; % Set this position accordingly to the problem. 

N=sqrt(k^2-(2*pi*fx).^2-(2*pi*fy).^2); % Look in chapter 4 for reference. 

Nm=sqrt((n*k)^2-(2*pi*fx).^2-(2*pi*fy).^2); % Look in chapter 4 for reference. 

num=4*N.*Nm*exp(-i*N*L); % Numerator of the equation. 

den=(Nm-i-N).^2.*exp(-i*Nm*L)-(Nm-N).^2.*exp(i*Nm*L); % Denominator of the 

equation. 

% Required terms for taking Fourier transform 
A0=exp(-i*2*pi*fx*x0).*exp(-i*2*pi*fy*y0); 


% Synthetic Data Generation 

% Transmittance for a plane wave through a given media 
t=4*n*exp(-i*k*L)./((n-i-1 )^2*exp(-i*n*k*L)-(n-1 )'^2*exp(i*n*k*L)); 


% scattering Points in the spatial frequency domain. 
xl=-l; yl=-l; zl=-l; 
x2=-l; y2=l; z2=-3; 
x3=0; y3=0; z3=-5; 
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x4=l; y4=l; z4=-1.3; 
x5=l; y5=-l; z5=-2.2; 


% Incident field after the penetration of the media. 
El=t*A*exp(-i*k*zl); 

E2=t* A*exp(-i*k*z2); 

E3=t* A*exp(-i*k*z3); 

E4=t* A*exp(-i*k*z4); 

E5=t*A*exp(-i*k*z5); 


% Position of the scatter points in frequency domain for the given 
% measurement plane at Z=constant 

% Since the scatter points are in terms of dirac delta functions, 

% the field measurement would become the summation of the Green's function 

% and the incident wave after performing the Rytov Approximation 

G1 =(num./den). *exp(i*N*Z). *exp(i*N*z 1) ./(2*N); 

G2=(num./den).*exp(i*N*Z).*exp(i*N*z2)./(2*N); 

G3=(num./den).*exp(i*N*Z).*exp(i*N*z3)./(2*N); 

G4=(num./den).*exp(i*N*Z).*exp(i*N*z4)./(2*N); 

G5=(num./den).*exp(i*N*Z).*exp(i*N*z5)./(2*N); 

% Required terms for taking Eourier transform 

A1 =exp(-i*2*pi*fx*x 1). *exp(-i*2*pi*fy *y 1); 

A2=exp(-i*2*pi*fx*x2).*exp(-i*2*pi*fy*y2); 

A3=exp(-i*2*pi*fx*x3).*exp(-i*2*pi*fy*y3); 

A4=exp(-i*2*pi*fx*x4).*exp(-i*2*pi*fy*y4); 

A5=exp(-i*2*pi*fx*x5).*exp(-i*2*pi*fy*y5); 

G=G 1 * A1+G2* A2+G3 * A3+G4* A4+G5* A5; 

% Generating Data 

D=E1. *(exp(G 1. * A1)-1 )+E2.*(exp(G2.* A2)-1 )+E3. *(exp(G3. * A3)- 
1 )+E4.*(exp(G4.* A4)-1 )+E5. *(exp(G5. * A5)-1); 

% Setting matrix dimension 
[MM,NN]=size(D); 

D=D(1:MM,1:NN); % Change matrix dimension if warranted. 

% Converting bin numbers to real axis numbers in spatial domain(meters). 

% Range has to be up to Nyquist frequency. 
ax=linspace(-fn/2,fn/2,length(D)); 

no=(randn(MM,NN))+i*(randn(MM,NN)); % noise in complex values in the frequency 
domain 

% set the signal to noise in dB 
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SNR=-20; 

S=abs(D.*conj(D)); 

NO=abs(no.*conj(no)); 

signal=sum(sum(S)); 

noise=sum(sum(NO)); 

factor=sqrt((signal/noise)*10^(-SNR/10)); 

no=no*factor; 

d=D+no; % data with noise in the frequency domain 
% End of (noisy) Data Generation 

% Tikhonov Regularization 
[u,s,v]=svd(G); 

for alpha=linspace( 10^-2,10^2,1000) % Range of Regularization parameters 

Reg=inv(alpha*eye(size(G'*G))+G'*G)*G'*d; % Regularized solution 
sol=norm(eye(size(G'*G))*Reg); % Solution norm 

res=norm(G'*G*Reg-G'*d); % Residual norm 


figure(3) 

loglog(res,sol,'-o'); 

hold on 

title('L-curve for Tikhonov Regularization'); 

xlabel('Residual norm'); 

ylabel('Solution norm'); 

figure(4) 

plot( alpha,res); 

hold on 

title('alpha VS residual norm'); 
xlabel('alpha'); 
ylabel('Residual norm'); 

end 

title('L-curve for Tikhonov Regularization') 
xlabel('Residual norm'); 
ylabel('Solution norm'); 
axis([10^-4 lOM 10^0 10^8]) 
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